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Abstract 

A recent Monte Carlo simulation determined the potential of mean 
force between two lysozyme molecules in various aqueous solutions [M. 
Lundetal. Phys. Rev. Lett. 100, 258105 (2008)]. The study involved 
a combination of explicit solvent and continuum model simulations 
and showed that there are significant ion-specific protein-protein in- 
teractions due to hydrophobic patches on the protein surfaces. In this 
paper we use the results of their study to determine the phase di- 
agram for lysozyme for aqueous solutions of NaCl and Nal. Two of 
the three phase diagrams have a stable fluid-fluid critical point, while 
the third has a slightly metastable critical point. This results from 
a secondary extremum in the potential associated with a repulsive 
interaction. This repulsive interaction reduces the effective range of 
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the attractive interaction and produces a metastable critical point. 
We compare the results of one of these phase diagrams with that for a 
model that includes ion-dispersion forces, but does not contain solvent 
structural effects. 



2 



1 Introduction 



Hofmeister effects refer to the relative effectiveness of anions or cations on a 
wide range of phenomena and date back to the original work of Hofmeister 
in 1888 [Tj. Examples include the surface tension of electrolytic solutions, 
bubble-bubble interactions, micelle and microemulsion structure and wet- 
tability p]. Hofmeister's original observation was that protein salting out 
depends on salt type, as well as on the salt concentration. In general, the 
effective interactions between most charged and neutral objects in aqueous 
solutions depend not only on the salt concentration, but also on the salt type. 
This fact has remained a challenge to theorists; a wide range of explanations 
have been proposed [3l H], most of which are not of a quantitative nature. 
Recently, however, there have been two parallel developments that have pro- 
vided considerable insights as to the causes of the Hofmeister effect. The 
first approach emphasizes the importance of including nonelectrostatic, ion- 
dispersion forces [21 El El [6] that had previously been neglected in standard 
treatments such as the Derjaguin-Landau-Verwey-Overbeek theory. In this 
treatment water is treated as a continuum, whose properties are described 
by a bulk dielectric constant and ion-dispersion forces between ions and so- 
lute particles are included, together with the standard Coulomb and van der 
Waals interactions to obtain a more self-consistent theory. In the second 
approach molecular dynamic simulations have been carried out to obtain 
the effective interactions between ions and interfaces, including the air-water 



3 



interface [3 [S] and the hydrophobic sohd- water interface P [ID]- In these 
MD simulations, water, surface and ion interactions are described by a phe- 
nomenological model (parametrized by Lennard- Jones interaction ranges and 
depths, and partial charges) whose parameters are chosen to match experi- 
mental bulk data. These studies have been successful in explaining a variety 
of ion-specific results, including explaining the adsorption of weakly hydrated 
ions such as bromide or iodide at the air- water interface [TJ [8] and the fact that 
large halide ions are attracted to hydrophobic solid surfaces, while smaller 
anions are repelled[TU]. These two complementary approaches have demon- 
strated the importance of including short-range interactions that account for 
both ion dispersion forces and short-range ion hydration effects. Indeed, a 
recent paper uses a parametrized potential of mean force that includes both 
ionic dispersion forces and short-range ion hydration to study ion specific 
effects for the double layer pressure between two uncharged interfaces 

The most recent application of molecular dynamics simulations has been 
to provide a step toward including a solvent induced ion-specific surface in- 
teraction in a Monte Carlo (MC) study of the interaction between lysozyme 
molecules in an aqueous salt solution[T2]. In this study the solvent is treated 
as a dielectric continuum, but solvent structure effects are implicitly included. 
The authors obtained the potential of mean force between two lysozyme 
molecules in various solutions of sodium chloride and sodium iodide, respec- 
tively. The corresponding second virial coefficients were shown to be in rea- 
sonable agreement with experimental data. They find that there are at least 
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two effects responsible for the Hofmeister series in this approach. Namely, 
it is the direct interaction of hydrated anions with positively charged amino 
acid residues as well as the affinity of these anions for hydrophobic patches 
at protein surfaces. The former interactions are stronger for chloride than for 
iodide, whereas the opposite is true for the second effect. Thus the effective 
protein-protein interaction in a particular salt solution results from a subtle 
balance between these (and perhaps other) effects. The fact that there is a 
stronger lysozyme-lysozyme interaction in aqueous Nal than NaCl is due to 
the hydrophobic effect of iodide surpassing the ion-pairing effect of chloride. 

An important property of an electrolytic solution of proteins is its phase 
diagram. It has been shown for a significant number of proteins that optimal 
crystal nucleation tends to occur when the solution is prepared in the vicin- 
ity of a metastable critical point of a protein-poor, protein rich fiuid-fiuid 
phase separation curve [T3l fl^ ITS] . Thus it is of interest to know the phase 
diagrams of the model of lysozyme obtained in the above molecular dynamic 
simulations. In this paper we present MC simulations of these phase dia- 
grams for several of the electrolyte solutions considered in reference [12] and 
compare the results of one phase diagram (0.2 M NaCl) with a corresponding 
calculation using a model based on ion-dispersion forces[5l [16]. We present 
in Section 2 a summary of the model and our simulation details. We give 
our results in Section 3 and present a brief conclusion in section 4. 
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2 Model and Methods 

The model of lysozyme in aqueous salt solution model, as given by Lund et 
al. [12] , consists of the following: each amino acid on the protein is represented 
as a soft charge-neutral sphere, located at the residue center of mass given by 
X-ray structure ALZT in the protein data bank. In addition, the protonation 
sites of all titratable groups are included at their original positions and their 
corresponding electrostatic charges are set for a pH of 4.7. The rigid protein 
molecules are allowed to rotate and translate in their MC simulations, while 
mobile salt and counter ions are explicitly included as soft spheres. The 
solvent is treated as a dielectric continuum and solvent structural effects 
(hydrophobicity) are implicitly accounted for. Specifically, the interaction 
energy between the ions and the nearest hydrophobic residues, Vspivij), is 
based on an empirical potential, from results obtained from MD simulations 
of ions in the presence of an air/water interface [T0l|T7]. The expression for 
the total interaction energy, from which the angular averaged potential of 
mean force (PMF) was calculated, is given as 




where Ib = 7.1 A is the Bjerrum length, rij is the distance between particles i 
and j, z is their charge number and e^j and the Lennard- Jones param- 

eters. From the above expression for PU, Lund et al. used MC simulations to 
calculate the angular averaged PMF, W{r), between the lysozyme molecules 
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for various aqueous electrolyte solutions. With their PMF data, we fit these 
potentials for the following systems: 0.2M NaCl, 0.3M NaCl and 0.2M Nal, 
as shown in Fig. [TJ It is important to note that the simulations carried out in 
Lund et al. were conducted at T = 298K = Tq. We use these potentials for 
the range of temperatures considered in this paper. For each model, we have 
defined the corresponding value of a to be the smallest value of r at which 
the potential crosses W = 0. For the three models examined ( 0.2M NaCl, 
0.3M NaCl and 0.2M Nal ), these values are nearly identical and are equal 
to 28.84A, 28.84A and 28.88A respectively. Therefore, for the remainder of 
the paper, a is to be interpreted as the a for the corresponding model being 
referenced. In addition, the potentials have been set equal to for values of 
r/a greater than 1.8, 2.06 and 2.02, respectively. 

2.1 Methods 

Using these PMFs, we perform MC simulations in order to determine the cor- 
responding phase diagrams and pair correlation functions. Our systems con- 
sist of N=500 particles in cubic simulation cells subject to periodic boundary 
conditions. The same number of MC cycles are performed for both equili- 
bration and production, although the total number varies depending on the 
type of simulation. A single MC cycle is defined as N=500 MC steps where 
a step is a random choice from the usual repertoire of MC moves. 
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2.1.1 Fluid-Fluid coexistence 

We use the Gibbs Ensemble MC method [HI [19] to obtain the equilibrium 
coexisting densities of the protein-poor and protein-rich fluid phases. This 
method avoids problems associated with the formation of an interface be- 
tween the dense and dilute fluid phases that would otherwise be present in 
single cell simulations. In this ensemble, two physically separate simulation 
cells are coupled to the same heat bath and are used to emulate the two fluid 
phases that are in contact. Standard particle displacements are performed 
within each simulation cell; in addition, volume and particle exchanges are 
performed between the two cells. These exchanges are chosen such that the 
total volume and number of particles of the system are conserved and the 
simulations obey detailed balance. On average, we chose the ratio of par- 
ticle displacements to volume moves to be 100:1; the frequency of particle 
transfers is chosen to give reasonable acceptance rates of approximately 1- 
3%. The equilibrium and production run times are at least 2* 10^ MC cycles 
each, where a MC cycle in this ensemble is N=500 attempts at one of three 
randomly selected trial moves: particle displacements, volume exchanges or 
particle exchanges. 

2.1.2 Fluid-Solid coexistence 

Fluid-solid coexistence are obtained via the Gibbs-Duhem method. This 
method involves integrating the flrst-order Clausius-Clapeyron equation ^ = 
'where P is the system pressure, P = 1/kT, A indicates a difference 
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between the liquid and solid phases and h and v are the molar enthalpies 
and volumes respectively. One caveat to this approach is that it requires 
the knowledge of an initial coexistence point on the /5P — jSfi plane. Con- 
sequently, we carry out a series of NPT simulations along an isotherm to 
determine the equation of state for both the fluid and solid phases. The 
equilibrium and production times for each NPT simulation are taken to be 
equal and at least 2 * 10^ MC cycles. Each isotherm requires a minimum of 
50 data points to obtain accurate fits. Once the equations of state are ob- 
tained, they are fit to the following form for the liquid pPLig = ^iiiZ^Y 
and jSPsoi = J2t=o ^iP" solid, where the a, Cj and m are distinct for 

each isotherm examined and chosen to best fit the data (Whereas in other 
studies, one typically use m=3, we find it necessary to increase it to 5 or 6 
in cases studied here.). Integration of the equations of state will yield the 
free energies and chemical potentials for the liquid and solid. However, as 
a result of the integration, one must know the free energy of the solid at a 
particular reference density along the solid isotherm. To calculate this free 
energy, we use the Frenkel-Ladd method for soft-core continuous potentials 
[20] to harmonically couple the particles to lattice sites. With the free en- 
ergy of the solid at the reference state known, straightforward integration 
produces the expressions for the corresponding chemical potentials and thus 
the initial coexistence points are determined via the conditions for mechani- 
cal and chemical coexistence in equilibrium, i.e. PuqipLiq) = Psoi{psoi) and 
PLiqiPLig) = Psoi{psoi) where pLiq and psoi are the coexisting hquid and sohd 
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densities. 



2.1.3 Accurate determination of Tc 

It is known that systems with short-range attractive interactions (< 1.25cr) 
give rise to hquid-hquid phases separations that are metastable with respect 
to the solubihty curve [2T] . i.e., freezing preempts fluid-fluid phase separa- 
tion. Choosing as an initial state a position near a metastable fluid-fluid 
critical point is desirable for experimentalists trying to grow high quality 
protein crystals from solutions as nucleation rates reach a maximum in that 
region |15]. An important observed characteristic of metstable fluids is their 
small and negative second virial coeflicients B2 flrst noted by George and 
Wilson[T3]. As will be discussed in more detail in our results, the phase 
diagram for the 0.3M NaCl system is only very slightly metastable as a con- 
sequence of an additional repulsive maximum in the potential. This is impor- 
tant because the system would otherwise exhibit a stable fluid-fluid transition 
without the addition of the repulsive maximum. Therefore in an effort to sub- 
stantiate that this system is indeed metastable, an accurate estimate of the 
critical temperature is determined by the Bruce- Wilding flnite-size scaling 
method [22]. The location of the critical point can be obtained by matching 
the probability distribution of the ordering operator M = ^5^7 with the uni- 
versal distribution characterizing the Ising class. The number density and 
energy density are deflned by p = L~'^N and u = L~''^U, respectively, with U 
the total energy of the system, d the dimension of the system. To obtain the 
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value of Tc{L) for a given simulation cell length L, we run a series of Grand 
Canonical Ensemble simulations and "tune" the values of chemical potential, 
temperature and s until the probability distribution of the ordering operator 
P{M) collapses onto the universal P*{M) for the Ising model. When good 
estimates of the tuning parameters are obtained, we perform longer runs 
of 2 * 10^ MC cycles to collect accurate statistics followed by a histogram 
reweighting procedure to obtain accurate values of the tuning parameters at 
criticality. This procedure is repeated for increasing values of L and then 
extrapolated to the limit L = oo to determine the critical temperature of the 
infinite system. A detailed discussion on this procedure can be found in Li 
et al. [23] . 

3 Results 

We have determined the complete phase diagrams for three different aqueous 
solutions of lysozyme, 0.2 M and 0.3 M NaCl and 0.2 M Nal shown in Figures 
His] and [6] respectively. In the case of 0.2 M NaCl, we obtain a "normal" phase 
diagram, in that the fiuid-fiuid coexistence curve is stable. In the case of 0.3 
M NaCl, we find that there is a metastable fiuid-fiuid coexistence curve, 
although the metastability gap (defined as (T^ — Tc)/Tc, where Tc denotes 
the critical temperature for the finite system, and Tl denotes the temperature 
on the solid curve that corresponds to the critical density) is very small. We 
estimate the triple point for the 0.2M solution to be approximately 0.4 Tq. 



11 



The situation for the 0.3M NaCl solution is a httle trickier, since the sohd 
branch of the sohd-hquid coexistence curve almost touches the fluid-fluid 
curve. In order to determine the critical point for the infinite system, we 
employ the method of finite size scaling (FSS), discussed briefly in section 
2.1.3. Our results for Tc{L) are shown in Fig. [2j (See the figure caption for 
more details.) The extrapolation L — > oo gives us a reliable estimate of the 
critical point, Tc(oo), for the bulk system. However, the NPT simulations we 
use to determine the solid-liquid curve are also subject to finite size effects 
and we are unable to extrapolate these to obtain their bulk limits. Therefore 
it is difficult to be precise about the magnitude of the metastability gap. 
Thus it would seem prudent to conclude that fluid-fluid critical point for 
0.3M NaCl is right on the edge of being metastable. In contrast to our 
results for this model, another model that takes into account ion-dispersion 
forces [IB] leads to a large metastability gap of 8.1%. 

For the 0.2M NaCl and 0.2M Nal solutions, we obtain estimates of the 
critical points by fitting the fluid-fluid coexistence data to the following stan- 
dard equations 

^^^p, + A|r-r,| (2) 

Pi-Pg-B\T-T,\^ (3) 

where Tc and pc are the critical temperature and density, respectively, and 
P ^ 0.326 is the 3D-Ising critical exponent. Finite size effects would have to 
be taken into account to obtain the critical point parameters for the infinite 
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system. In the case of 0.2M Nal ( Fig. 0), we estimate the triple point tem- 
perature to be around 0.7 Tq. We are not able to make direct comparisons 
between simulation results and experimental data for either salt, however, 
since to the best of our knowledge, no experimental data exist for the mo- 
larities studied in this paper. However, for solutions with 0.5M or greater 
concentrations of NaCl, experimental data is available. This data shows the 
presence of a metastable fluid-fluid phase for all salt concentrations greater 
than 0.5M [2l]. It is also likely that there is a metastable critical point for 
concentrations somewhat smaller than 0.5M, given that the NaCl solution at 
that molarity exhibits a large metastability gap |24] . 

It is also interesting to note that for both stable phase diagrams ( Figs. 
Il]and[6]), there is a metastable continuation of the liquidus curve along the 
high density branch of the fluid- fluid coexistence curve below the triple point, 
Ttp. Namely, as we continue our simulation of the coexisting liquid phase in 
the liquidus line below Ttp, this becomes the metastable dense liquid phase 
we have obtained by our gibbs ensemble method. 

Figure [3] plots the pair correlation functions of the liquid and solid at 
conditions corresponding to phase coexistence for the 0.2M Nal and 0.3M 
NaCl lysozyme solutions, respectively. In order to show that the solid does 
indeed exhibit long range order, we ran simulations of N=4000 particles to 
probe larger distances. It is clear that both liquid and solid correlation 
functions display an oscillatory behavior. The reason that it does so for the 
liquid phase is due to its high density. Note that in Figures H] and [6] the 
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separation between the liquid and solid coexistence densities is very small. 
However, the correlation functions for the solid still shows a strong peak at 
the largest distances, suggesting a crystal structure, while it already decays 
to 1 for the liquid phase after a few (~ 4.5) as. 

One can see from the phase diagrams of Figs. [5] and [6] that the details of 
the corresponding potentials ( Figs. [T](b) and [U^c) ) are important in deter- 
mining the phase diagram. Although the attractive wells of their potentials 
are similar in both the width [25] and the depth, their phase diagrams are 
very much different. The 0.3M NaCl solution has a metastable fluid-fluid 
coexistence curve while the 0.2M Nal solution has a stable phase separation 
curve. This can be largely attributed to the fact that the potential for the 
0.3M NaCl has a repulsive region which reduces the effective range of the 
attractive interaction. In other words, the critical point of the fiuid-fiuid 
coexistence curve is driven to a lower temperature by the additional positive 
extremum in Fig. [I](b). This behavior was also observed by Brandon et al. 
[26] . in their study of models with multiple extrema. They argued that for 
more complicated potentials than those with a single extremum, the effec- 
tive range of attraction defined by Noro and Frenkel[27j was more useful in 
characterizing the phase diagrams than, say, the width of the attractive well. 
Noro and Frenkel[27] developed an extended law of corresponding states, in 
which they introduced an energy scale (the depth of the attractive well), an 
effective hard core cxe// and a range R. To do this for systems whose poten- 
tials are continuous, one divides the potential into attractive and repulsive 
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terms, Vatt and Vrep, respectively, and defines an "equivalent" hard-core di- 
ameter for the repulsive part of the potential using an expression suggested 
by Barker and Henderson |28j: 



^0 

The reduced second virial coefficient is then defined as B2 = -B2/(— -y^)- 
One then defines an effective range R of the attractive interaction by equat- 
ing B2 of a square well system with that of the system in question. (Note 
that the range of the attractive interaction for a square well system is defined 
unambiguously, whereas this is not the case for other models.) Although this 
range is temperature dependent, it provides a useful length scale for models 
such as those that have been used in colloidal and protein solutions. The 
value of this reduced second virial coefficient at the critical point has been 
found to be relatively constant for a wide range of models that are com- 
monly used to describe phase transitions, ranging from extremely narrow 
attractive wells (Baxter's adhesive hard sphere model) to the van der Waals 
limit of infinitely long-range attractive wells[27]. These values are in the 
range —2.36 < i?2(7"c) < —1.33. We show the reduced second virial co- 
efficient B2{Tc), the effective hard core diameter a^jf and the range R for 
all three potentials studied here in Table [H All our values for B2 (Tc) are 
in the same range as for the models cited in Ref. [27]. Noro and Frenkel 
also estimated that the fiuid-fiuid transition became metastable with respect 




(4) 
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to freezing when R = 0.14, consistent with a variety of models that had 
been studied at that time. Subsequently, however, it was shown that for the 
square well model R ^ 0.25[21J. This latter value has been proposed based 
on a simple van der Waals model for both the fluid and the solid phase [29] 
and by a phenomenological argument using a cell model for the solid [30]. If 
we include the region of repulsive interaction associated with the secondary 
maximum in our calculation of the effective hard core diameter, we find that 
the range for the 0.3 M NaCl system that is slightly metastable is i? ~ 0.20. 

Finally, we note that we have also run preliminary NPT simulations for 
the 0.3M NaT solution. However, no phase diagram was calculated due to 
the fact that the difference in density between the liquid and solid isotherms 
is extremely small even when the pressure is sufficiently high. In addition, 
at higher pressures, the density on the liquid isotherm has large fluctuations 
about its average; in fact, it typically transitions from the liquid phase to 
the solid phase where it remains stuck for the remainder of the simulation. 
Therefore, since we are unable to distinguish sufficiently accurately between 
the liquid and solid phases, we are unable to perform an accurate calculation 
of the 0.3M Nal phase diagram by this method. 

4 Conclusion 

We have studied models of three electrolyte solutions of lysozyme and demon- 
strated the effects of hydrophobic surface forces on the phase diagram of 
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lysozyme. Although the molarities studied are too small to see a significant 
metastability gap, the 0.3 M NaCl solution does have a slightly metastable 
fluid-fluid critical point. Interestingly, the small range R associated with 
this is produced by the effect of a repulsive interaction associated with a sec- 
ondary extremum. In the absence of that repulsive region the system would 
have a stable fiuid-fiuid transition. This is another example of a phenomenon 
first studied by Brandon et al. j26j. It is not clear what the physical origin 
of the secondary maximum is for the model of Lund et al. [12], but Brandon 
et al. PB] attribute their secondary maximum to the effect of water restruc- 
turing near the solute particles. It would be interesting to have a better 
understanding of its origin for the current model. As we are unable to com- 
pare the results of this study with any experimental results for lysozyme, we 
cannot determine the accuracy of the model. It would seem likely that fur- 
ther developments will include in addition the specific ion-dispersion effects 
discussed by Bostrom et al. [5J. These will require more accurate quantum 
mechanical calculations of the amplitudes of the ion- dispersion forces; such 
calculations are currently being carried out [311 132] . 
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Figure 1: Total potentials for the (a) 0.2M NaCl (b) 0.3M NaCl (c) 0.2M 
Nal electrolyte systems. Solid line is the fitted polynomial and dots are data 
obtained from MC simulations. 
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Figure 2: FSS results for Tc(L) as a function of L-^^+^^^Z" for the 0.3M NaCl 
electrolyte solution, v is the critical exponent for the correlation length and 9 
is the universal correction to scaling exponent. For the 3D Ising universality 
class, V = 0.629 and (3 = 0.326. By extrapolating to infinite volume (L — >■ 
oo), we can obtain an estimate of the true bulk behavior: Tc(oo) ^ 0.755To. 
For more details of FSS, see Ref. [23] 
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Figure 3: Plot of the liquid and solid pair correlation functions corresponding 
to coexisting liquid and solid states for (a) 0.2M Nal at T = 1.6To,P = 
20A09(r^/kTo,pLiq(r'^ = 1.069, psoia^ = 1-104 and (b) 0.3M NaCl at T = 
0.91To,P = 2m2aykTo,pLiga' = 0.75A, psoia^ = 0.798. 
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Figure 5: Phase diagram curve obtained from Monte Carlo simulations for 
the aqueous lysozyme solution with NaCl electrolyte at 0.3M. The fluid-fluid 
separation curve is metstable. Open circles are the liquid-solid coexistence 
data and solid diamonds are the fiuid-fiuid coexistence data. The solid star 
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